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Abstract. 

The possibility of using hexagonal structures in general and graphene in particular 
to emulate the Dirac equation is the basis of our considerations. We show that Dirac 
oscillators with or without restmass can be emulated by distorting a tight binding 
model on a hexagonal structure. In a quest to make a toy model for such relativistic 
equations we first show that a hexagonal lattice of attractive potential wells would 
be a good candidate. First we consider the corresponding one-dimensional model 
giving rise to a one-dimensional Dirac oscillator, and then construct explicitly the 
deformations needed in the two-dimensional case. Finally we discuss, how such a 
model can be implemented as an electromagnetic billiard using arrays of dielectric 
resonators between two conducting plates that ensure evanescent modes outside the 
resonators for transversal electric modes, and describe an appropriate experimental 
setup. 
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1. Introduction 

Dirac operators play a central role in relativistic quantum dynamics, from the early 
work of Dirac and the exact solutions for the hydrogen atom to later work including 
quantum chromodynamics and the Dirac oscillator [H Ej |3]. Based on early work of 
Wallace there has been much recent work centered around the fact that the mean 
field theory of grapene is well represented by the free Dirac equation near the edges of 
the first Brillouin zone, i. e. at the center of the band with so-called Dirac points around 
which linear dispersion relations hold O [6l [7]. Simulations of this kind of situation 
are ongoing at various labs using hexagonal and occasionally triangular arrays in single 
particle problems or classical waves, including acoustics [8j, microwaves and photonic 
crystals [9] as well as true quantum simulations in nanostructures [TOl [TTl [121 [IS]- We 
wish to recall, that any lattice with coordination number three will yield points with 
approximate linear dispersion relations similar to Dirac points, but we need a tight 
bindig situation to guarantee isotropic cones and a hexagonal structure to introduce an 
additional discrete degree of freedom for the small and large component of the effective 
spinor. This becomes very transparent as we note, that using a hexagonal structure 
made up of two different triangular structures, we obtain a finite gap in the spectrum, 
i.e a Dirac equation for a massive particle, a situation that would correspond to a B-N 
(Boron Nitride) lattice. 

In this paper we will adress the task, to find similar analoga for Dirac operators, 
that describe situations other than the free particle. Formally we shall assume arrays of 
potential wells, that can hold exactly one bound state. In between the wells these states 
decay exponentially. As the overlap of functions of two wells describes the coupling, 
this implies to a good approximation a tight binding system. In this framework in 
principle one and two - dimensional Dirac type problems can be formulated, and in 
some cases exactly solved. We shall here focus on the Dirac oscillator as introduced 
by Marcos Moshinsky; this system is frequently used [HI [151 [HI El [HI IHl [20] , and is 
particularly simple to implement. Other solvable problems, such as gyroscopes |21] and 
the Coulomb problem as well as systems with random potentials will be touched upon 
in the conclusions. 

In particular, the one and two dimensional Dirac oscillators have been considered 
recently in the context of both relativistic quantum mechanics and quantum optics 
[22l [231 El]- Their importance as paradigmatic integrable models is well established 
and their realization in the context of tight binding models is desirable. We want to 
stress that the dimensionality in our examples plays a crucial role in their algebraic 
properties and spectra. As we shall see, the exceptional infinite degeneracy of the two 
dimensional Dirac oscillator (as opposed to the finite degeneracy in one dimension) 
is intimately related to its specific realization on the lattices described above. It is 
therefore interesting to consider both dimensionalities, as they exhibit clear differences. 
It remains as an open question, whether three dimensional relativistic wave equations 
can be emulated on a three dimensional lattice. 
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The description of such systems will depend on a massive distortion of the hexagonal 
lattice. This causes probably prohibitive obstacles to any implementation in terms of 
mean fields of graphene related structures. For recent theoretical work dealing with 
deformations on carbon sheets see [25]. For classical wave models on the other hand 
such distortions can be implemented and we shall discuss specifically how this is done in 
what we usually call a microwave billiard. In two dimensions such systems yield scalar 
wave equations, and can readily be interpreted as single particle quantum systems [26] . 

Microwave billiards have been used to simulate a wide variety of phenomena 
including quantum chaos [271 ISH [291 [30] , scattering from open billiards [3T1[32], transport 
phenomena [331 El] , fidelity decay [351 EE] and disordered systems [33, [3H] , as well as 
recent work on Dirac points [39], [H]. This wide range of success combined with recent 
experiments with arrays of small dielectric micro wave resonators gives us good reason 
to hope for interesting results in orderd structures such as the ones we need. The 
evanescence of the cavity modes, that we intend to use, guarantees exponential decay 
and thus tight binding situations can be emulated. 

In the next section we study a one dimensional crystal; this will serve us to fix 
notation and basic ideas in a simpler context, that is not without relevance in itself. 
This toy model reveals Dirac-like hamiltonians in both periodic and deformed structures. 
In section 3 we analyze and reproduce similar results for the two dimensional case 
with and without deformations. We also consider corrections to nearest neighbour 
interactions and obtain the form of energy surfaces beyond tight binding. In section 4 
we make contact with experimental applications of our treatment by considering arrays 
of resonators in microwave cavities. The applicability of the tight binding model in 
this context is discussed. Finally, we draw some conclusions and give an outlook on a 
selection of other Dirac systems, that can be emulated. Useful results are included in 
the appendix. 

2. One-dimensional crystal 

In this section we start out with two lattices. First, we consider the Schroedinger 
equation with square wells supporting a single bound state, located in a periodic one 
dimensional array. Here, traditional aspects of the existing theory {i.e. Bloch waves) are 
reviewed under an algebraic approach. The specific form of the localized wave functions 
is ignored, as it results to be irrelevant. An analogy between degeneracy points of the 
spectrum and a one dimensional Dirac point is revealed in this toy model. Near these 
points we shall cover a one- dimensional Dirac equation. 

Once this is achieved, we proceed to deform the lattice from its periodic 
configuration by forcing a specific operator algebra, namely the one corresponding to 
harmonic oscillator ladder operators. As an outcome, a one dimensional Dirac oscillator 
shall be realized. 
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Figure 1. Configuration of potential wells (or resonators) on a one-dimensional 
lattice. Chain a) shows the periodic case. Chain b) corresponds to a general 
deformation of both sublattices. Chain c) shows the resulting deformation after 
imposing the harmonic oscillator algebra. 



2.1. The Dirac point in one dimension 

We start by fixing some notation. In the following, we adopt natural units h = c = 1. 
A lattice consisting of two periodic sublattices is considered. Let A be the distance 
between neighbouring potential wells. They have the same period and are denoted as 
type A and type B. Each sublattice point can be labeled by an integer n according to 
its position on the line, i.e.Xn for type A and for type B (see figured]). The energy 
of the single level to be considered in the well is denoted by a and (3 for type A and 
B respectively. The state corresponding to a particle at site n of lattice A is denoted 
by \n)A and the corresponding localized wave function is given by ^a^x — Xn) = (x|n)^. 
For lattice B we define the wavefunction E,b{x — yn) = {x\n)B. Our hamiltonian is well 
approximated by a tight binding model if the overlap between localized wave functions 
can be neglected and the nearest neighbour coupling matrix element A is taken as the 
first off-diagonal element of the hamiltonian in the localized basis. 

For convenience we have split the lattice into two sublattices and we write the 
hamiltonian in the form 
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H 



Haa Hab 
Hba Hbb 



(1) 



rather than in the usual tridiagonal form. The entries of each block are given by 
H^j^ = {n\Hij\m) with n,m integers. Thus, each block extends from {—oo\Hij\ — oo) 
to {oo\Hij\oo) . In this basis we write the hamiltonian as 



/ . 
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a 



A A 
A A 
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A A 
A 



/3 
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(2) 



where the elements outside the indicated diagonals are all zero. Expression (j2]) can be 
cast in terms of Pauli matrices Us, cr+ = ai + ia2, cr_ = a\. by defining 11 = Hab ® 1 and 
setting M = {a- /3) /2, Eo = {a + /3) /2. We have 



H = Eo + asM + (T+n + (T_n^ 



(3) 



displaying explicitly the Dirac-like structure of the hamiltonian. To ensure the analogy 
between the Dirac hamiltonian and (|3]), we consider the spectrum of H. We note that 
n has the remarkable properties 



[n,tf] = o, nnt = A(n + nt). 

From these relations we may compute the spectrum by squaring H 



(4) 



{H - Eof = M^ + UUl (5) 
Note that Bloch's theorem manifests itself in the spectrum and eigenfunctions of Iin^ 



as 



U(j)k = A(l + e 



i2TTXk\ 



nnVfc = A^ii + e 



i2n\k\2. 



(6) 



with 
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n) = 



3i27rAfc(n+l) 



V 



(7) 



Therefore energies and eigenfunctions of H are given by 



(f>k 

Eo-M 



(8) 



where A?^ is a normahzation constant. When M = le. when the sites A and B are 
equal, we can find conical points of the form ko = ±1/(2A). Expanding around these 
points we find the usual expressions for the relativistic energies of Dirac particles with 
momentum k = k — ko, effective speed of light A: 



E(k) ^Eo± VA2«;2 + (9) 

This is vahd also for non-zero rest energy M in the case of different lattices. Then we 
find a gap in the spectrum. The eigenfunctions satisfy the Dirac equation in momentum 
space 



{E - Eo)iP^{k) = [aiK + a^M]^^{K). (10) 

In the next subsection wc shall proceed to show that we can obtain the Dirac oscillator 
by deforming the double lattice. 

2.2. One- dimensional Dirac oscillator and lattice deformations 

Using the notation of the previous section, we modify the energy operator by deforming 
the lattice of potential wells. This implies abandoning the periodic structure and leads 
to a site dependence of the couplings in the corresponding tight binding model. We 
denote by An^n+i the couphng between sites at yn and x^+i, while A^^n denotes the 
coupling between sites x„ and y„. These couplings are approximately proportional to 
the overlap between neighbouring sites. They decay exponentially as a function of the 
separation distance between the potential wells, i.e. 



A„,„+i = Ae-*-/^, A„,„ = Ae-</^ (11) 

where dn and d'^ are the deviations from the periodic configuration, i.e.dn + A = 
|?/„+i — a;„|, d'^ + \ = \yn — Xn\- When dn = d[^ = 0, the periodic configuration is recovered 
(see figure). The length A is the penetration depth into the classically forbidden region 
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for the wave function of the single welL We expect a modification of the operators 11, 11^ 
caused by the change in the position of the wells. One should keep in mind that such 
deformations have the effect of breaking the periodic symmetry of the system and Bloch 
wave functions cease to be useful. However, if one finds an exactly solvable model for 
a deformation which is continuously connected to the periodic configuration, then the 
corresponding solutions could be considered to constitute a generalization of Bloch's 
waves. 

The algebraic properties of observables in hamiltonians have a clear connection 
with integrability and exact solvability. Therefore, the simplest way to extend our 
hamiltonian is by replacing the operators 11,11''^ by a,a\ such that their algebraic 
properties are those of known solvable systems. In particular, we propose the harmonic 
oscillator algebra, since it is a paradigmatic example |43]. The hamiltonian (IHl) becomes 



H = Eq + a^M + cr+a + cr_a 



t 



(12) 



We require that such extensions reduce to the periodic case in some free limit. We 
propose then 



V 

and impose the condition 



^n+l,n+l ^n+l,n+2 



(13) 



[a, a^] = wA = constant (14) 

where u) stands for a frequency. If a; A = we recover the algebra of 11, 11^ (Bloch limit). 
Computing the commutator in ( |T^ . one finds the conditions 

where the second equality is a recurrence equation. The corresponding solution is 



and thus 



A 



n,n+l 



a/A2 - nuA, neZ 



(16) 



A VA2 - nuA 

A JA2 - (n - l)a;A 



V 



(17) 
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The first condition imphes c?^ = (a chain of dimers). The second condition relates the 
displacement between resonators dn with the corresponding couplings. 
For convenience we impose a third condition, namely the inequality 



dr, > d. 



(19) 



This condition ensures that the displacement dn increases monotonically as a function of 
the distance from an arbitrarily chosen origin. Combining it with the second condition 
this implies a finite grid, which in turn leads to a finite Hilbert space of dimension 
^max + 1 with Umax = [\~\]- The Operator a takes the finite matrix form 



f A v/A^^ 
A 

V 



7,.tWA 



A2 - (n. 



DujA 



(20) 



In principle other non-monotonical choices of scaling dn are possible and lead to other 
finite or infinite arrays. For our modeling purposes finite arrays are advantageous, and 
in any case our hamiltonian will emulate a Dirac oscillator only in the vicinity of zero 
energy (energies near Eq) as we will see later. 

We now have to check the commutation relations between a' and a'^ 



/ 1 



[a', a'^] = wAl + [A^ - (n. 



uAl + 0(l/n„ 



, + l)c^A] 



(21) 



Note that the correction term is of order l/umax as the principal term is accompanied 
by the identity, while the correction acts only in the first component of this basis. The 
prefactor is of order 1 since A/u ~ [|A/a;|] + 1. The distortion is 



dn = A log 



A2 



A2 - nuA 



< n < nn 



(22) 



Now we can replace a' by a and return to the hamiltonian (fT2|) for the finite system. 
Our construction allows to compute the spectrum and eigenf unctions up to a shift 
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smaller than unity. As we shall see, they correspond to those of the Dirac oscillator 
[3]. We consider eigenvectors 0„ of the number operator such that a^acpn = {^'^)n(l)n, 
^Vn = \JujA{n + l)(f)n+i and a^o = 0. The hamiltonian has an integral of the motion 
given by the operator / = a^a + |u;A((T3 + 1) with eigenstates 0n+i|— )• In this 

basis, H — Ef) is reduced to 2 x 2 blocks of the form 

( ^ ^u:A{n + l) \ 

\ ^uA{n + l) -M ) ' 

These blocks can be diagonalized, leading to the energies 



E{n) =Eo± ywA(riTl)TM2, 0<n< n^^x- (24) 

An additional 1x1 block is due to the singlet t/^o = 0o|— ), leading to if^/'o = (-E'o — M)-?/'o- 
The eigenfunctions corresponding to the doublets are obtained in the form 



^^+i = ArU„|+)+ ^ y)^ }^ 0„+i|-) , ^<n<n^ax. (25) 

V ^luJA{n + 1) j 

where iV is a normalization constant. The specific form of 0^ is given in the appendix. 

In the next section we shall see that a similar construction occurs very naturally 
if we deform hexagonal lattices in two dimensions, which in the undeformed case are 
closely related to the mean field theory of graphene and boron nitride sheets. 



3. Two-dimensional crystal 

The concepts given in the last section are now extended to two dimensions. We shall 
consider hexagonal structures that emulate two-dimensional Dirac equations and can 
emulate the mean field of graphene or Boron-Nitride near the gap at the center of the 
usual band. We shall use the same algebraic strategy to derive spectra and a possible 
extension through deformations, namely the two-dimensional Dirac oscillator. 



3.1. Hexagonal lattice 

Let us fix the notation for this system. The honeycomb lattice is divided in two 
triangular sublattices, one of them generated by the set of vectors ai = (3v^, 0),a2 = 
(-3/^2, 3/2), ag = (3/^2,-3/2) (type A) while the other sublatt ice is obtained by 
adding the vectors bi = (0, 1), bs = (-3/v^, -1/2), bg = (3/^2, -1/2). These vectors 
are given in arbitrary units (see figure [2]). We denote the linear combinations of a^ 
by A, where A is a vector parametrizing the points of sublattice A. For sublattice 
B we use the vector parameter A -|- bi. The position vectors r^,r5 of the periodic 
lattices are obtained by introducing the factor A. For periodic arrays this means 
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Figure 2. Vectors in a two dimensional hexagonal lattice. The vectors bj 
connect sublattice A with sublattice B. The vectors connect points in the 
same sublattice. 

Fyi = A A, = AB. In further considerations this notation will be useful, since deformed 
lattices admit a parametrization by vectors A, B, but the corresponding position vectors 
tajTb become more complicated functions of aj,bj. The state vectors for individual 
potential wells on grid A shall be denoted by |A), giving wave functions of individual 
wells as ^a(i" — i"a) = For grid B we use |A + bi). As before, we consider different 

energies for the wells on grids A and B. 

Similarly as in the one- dimensional case, the hamiltonian in the tight binding 
approximation is constructed as 

i/ = Q;^|A)(A|+/3^|A + bi)(A + bi| (26) 

A A 

+ A(|A)(A + b,| + |A + b,)(A|) 

A,i=l,2,3 

where the first two terms indicate the total on-site energy on grids A and B respectively, 
while the last sum indicates the nearest neighbour interaction with coupling strength 
A. We shall analyze this system by considering again a subdivision of the Hilbert space 
according to sublattices A and B. Due to the coordination number in this lattice, the 
matrix representation in section 2 is no longer feasible. However, we may construct the 
usual Pauli operators through the definitions 

<7+ = ^|A)(A + bi|, a_ = ai (27) 

A 

<73 = E|A)(A|-|A + bi)(A + bi|, l = ^|A)(A| + |A + bi)(A + bi|, 

A A 
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while the operators 11, 11^ are defined through 

n = ^ A (I A) (A + b, - bil + I A + bi)(A + b,|) . (28) 

A,i 

They have the algebraic properties 

[(T3, (T±] = ±2a±, [n, nt] = 0, [n, a,] = o (29) 

With M and Eq given as in section 2, we obtain the hamiltonian 

H = Eo + Mas + (x+n + (T_n^ (30) 

and 

{H - Eaf = M^ + ntf (31) 
Using Bloch waves, we have eigenvectors of the form [12] 



|A), 0t = Ee"''"^'"^''+'''^|A + bi) (32) 
for grids A and B respectively. They satisfy 



A.,i A,i 



nnVf = A^l ^ e*^.^^' '^! (33) 



'k 



The spectrum and the eigenfunctions are then 



E{k) = Eo± /a2| ^e^2^^b,.k|2 + (34) 



The degeneracy points of the spectrum for the massless case are kg = ±2^(1,— a/3). 
Expanding around such points one finds 



E(k - ko) - ^0 = ±VA2p + M2 (36) 
as expected. One can verify that the Dirac equation is again satisfied. 
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3.2. The importance of the tight binding approximation 

We have formulated a theory describing the propagation of waves in hexagonal arrays 
of resonators with and without deformations. The treatment has been successful in 
predicting the existence of Dirac points, in compliance with the common knowledge 
about this system when periodic symmetry is present. However, one may ask whether 
the tight binding approximation is an essential ingredient for a Dirac-like behaviour of 
the propagating wave. 

To answer this, let us review some of the assumptions we made and the 
corresponding properties obtained for the two-dimensional Dirac wavefunction. The 
spin emerges as the probability of the electrons to be located at sites of the triangular 
sublattice A (spin up) or B (spin down) . The momentum (or wave vector) as a conserved 
quantity came directly from Bloch's theorem, i. e. the periodicity of the system. Linearity 
around degeneracy points (which is essential to produce a Dirac hamiltonian) is a 
consequence of the hexagonal structure, while the existence of such degeneracy points 
came from the symmetry under the interchange of the two sublattices. We have seen 
that the appcarence of mass corresponds to the lift of such degeneracy. 

In addition to all these properties, one should consider isotropy as a fundamental 
requirement for the emulation of a free Dirac particle. We claim that rotational 
symmetry around degeneracy points is a direct consequence of the tight binding 
approximation, as we shall see. It is well known that rotational symmetry in the Dirac 
equation demands a transformation of both orbital and spinorial degrees of freedom. It 
is in the orbital part that we shall concentrate by studying the energy surfaces around 
degeneracy points beyond the tight binding model. 

Let us recall that the transition amplitudes between nearest neighbours (hopping 
transition) gave rise to the hamiltonian 



where the kinetic operator 11 could be constructed in terms of translation operators 
between nearest sites, i.e. 



with Tb- the one-site translation operator in the direction of bj. Since degeneracy points 
are located using the condition HipQ = for a Dirac state ipo, and given the fact that 
[n, n^^] = 0, it suffices to impose II'^o = 0. The operators Tb. are unitary and commute 
with each other, implying that their eigenvalues can be found simultaneously as e*^'*''*, 
where k is any real wave vector. Small deviations from degeneracy points (denoted by 
ko) in the form k = ko -|- k give the energy 



H = cr+n + (7_n^ 



(37) 



n = A E 7b,, 



(38) 



j=l,2,3 




(39) 
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which is rotationally invariant in k. However, one may try to introduce interactions 
between sites separated by more than one step under lattice translations. It is clear 
that a second-neighbour interaction of strength A' modifies the kinetic operator 11 as 



n = A Y ^b^+A' Y (ra.+T_.J, (40) 

i=l,2,3 1=1,2,3 

where the vectors aj have now appeared, connecting a point with its six second neighbors. 
The energy equation becomes 



E= |A^exp(zAk- bi) + A'^2cos(Ak- ai)|. (41) 

i i 

We expect a deviation of degeneracy points kj,, for which k = kg + K. Upon linearization 
of the exponentials in k we find the energy 



^(k-u)2 + (k- v)2 (42) 
where the vectors are given by 



u = AA ^ cos(Ak[) • hi)hi (43) 

i 

V = AA ^ sin(Ak[, ■ hi)hi + 2AA' ^ sm{Xk'f^ ■ ai)ai (44) 

i i 

Thus, the presence of A' yields the energy surfaces (H2l) as cones with elliptic sections 
whenever k, is inside the first Brillouin zone. Regardless of how we complete the energy 
contours to recover periodicity, it is evident that the resulting surfaces are not invariant 
under rotations around degeneracy points. The circular case is recovered only when 
A' = 0, leading to kg = ko. In this case, the vectors reduce to v = (1,0), u = (0, 1) 
when ko is the degeneracy point at (1/2A,0). 

In summary, extending the interactions to second neighbors has the effect of 
breaking the isotropy of space around degeneracy points, which is an essential property 
of the free Dirac theory. 



3.3. Two-dimensional Dirac oscillator 

Before analyzing lattice deformations, let us recall [191 122] that the two dimensional 
Dirac oscillator hamiltonian is quite similar to the one dimensional case, except for the 
replacement a h-j- a^, where = ax + iay is the chiral (right) anhilation operator in 
terms of cartesian anhilation operators a^^ay. The corresponding number operator is a 
conserved quantity and is given by Nr = N — L, i.e. the difference between the total 
number operator and the orbital angular momentum. Since the chiral left operator 
is absent in the expresion, the spectrum is infinitely degenerate. Eigenfunctions are 
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constructed as a combination of two usual harmonic oscillator functions with defined 
orbital angular momentum. Our aim is to produce the spectrum for this problem in 
the hexagonal array, together with a deformation which allows localization of the wave 
functions around some center. 

We proceed to deform the lattice through an extension of the kinetic operators, in 
analogy to the one dimensional case. Let us consider site dependent couplings A (A, A + 
bj) connecting the sites labeled by A, A + bj. Again, these are related to distances 
d{A, A + bj) between potential wells as A(A, A + bj) = Aexp(-rf(A, A + bj)/A). We 
define the ladder operator 

aR = J2 A(A, A + bi) (|A)(A + b, - bi| + |A + bi)(A + bj|) (45) 

A,i 

and impose [a^, a|j] = uA. After some algebra, one can prove that this leads to the 
conditions 

A(A,A + bi) = A, (46) 

A2(A, A + ba) + A2(A + b2, a + b2 - bg) = (47) 
A2(A + bi, A + bi - bg) + A2(A + bi - bg, A + bi + bg - bg), 

A2(A,A + b2) + A2(A,A + b3) = (48) 
A2(A + bi, A + bi - ba) + A\A + bi, A + bi - ba) + uA. 

The vector bi in the first equation was chosen arbitrarily, but due to symmetry a choice 
of b2 or ba would be equivalent. Due to the coordination number three, we obtain 
three relations rather than the two of the one-dimensional case. As complicated as 
the recursion relations may seem, one can easily construct a lattice reproducing them 
consistently. The relations ( H6i) and ( H7I) establish an equality between the lengths of 
opposite sides of a given hexagon. The relation (148|) containing u gives the deformation 
and can be split in two parts 

A\A, A + b2) - A2(A + bi, A + bi - b2) = uA sin^ 6 (49) 

A\A, A + bs) - A2(A + bi, A + bi - bs) = uA cos^ 6 (50) 

where 6 is an arbitrary angle. Having chosen previously the privileged direction bi, the 
angle 6 determines the relative stretching between directions b2 and on the seminal 
cell. The choice 6 = produces deformations only in one direction of the lattice (ba). 
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resembling the one dimensional case discussed above. This leads to a logarithmic law for 
the deviation distance similar to f lTTj) . Logarithmic stretching will hold for an arbitrary 
angle 6. To construct the grid in the general case, one starts with a regular hexagon 
as a seed and completes the scheme in figure [3] by extending lines of equal length in 
the direction of bi. Then, one completes the hexagons by drawing parallel lines for 
the opposite sides as shown also in the figure [31 Hexagonal cells satisfy the recursion 
relations above trivially. 

A restriction to a finite dimensional space occurs in a way similar to the one 
dimensional case. The maximum number of levels is A^^^ax = [|^|]- 

Considering a^j as the chiral operator restricted to this finite dimensional space, the 
resulting hamiltonian of this problem is 

H = Eo + a-iM + a+aR + (x^a^. (51) 

As the hamiltonian f lSTj) is formally identical to that in the one-dimensional case, we 
find the eigenvalues 

E^{Nr + 1) = Eo± ^uA{Nn + 1) + M\ 0<Nn<A/cu, (52) 
E{0) = Eq-M. 

The shape of the eigenf unctions is obtained by solving Oij^o = and applying the raising 
operators similar to appendix A. 

The hamiltonian does not depend on the left operators 0^,0]^ where = {aR)*. 
In the full space this would imply a Landau electron-like infinite degeneracy. The 
degeneracy does not occur for a fixed array of resonators, but it can be interpreted to 
reflect the arbitray choice of 6 if we consider the hypothetical use of an ensemble of 
arrays for all angles 6. 

Note that there is a physical limitation to the allowed degree of distortion which 
results from the fact that the nearest neighbours can change. Before this happens, the 
coupling of potential wells can no longer be dominated by the original three nearest 
neighbours. 

In figures HI |5l [6] we give some realizations of lattice deformations following the 
procedure indicated previously. The examples are related to the choices of ^ = 
0,7r/4, 7r/2, which determine the vectors to be deformed with a logarithmic law near 
the X axis of the graphs. 

4. Experimental implementation in microwave cavities 

To emulate Dirac-like equations on two dimensional arrays of resonators, we have to 
generate a scalar field. For electromagnetic fields, this can be achieved in a 2D metallic 
cavity which supports two types of independent modes: transverse magnetic (TM) mode, 
ip{r) = Ez{r), and transverse electric (TE) mode, ip{r) = Hz{r), r lying in the plane of 
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Figure 3. Construction of a deformed lattice satisfying the oscillator 
constraints. Here we show the direction in which the length of hexagons grow 
as a logarithmic function, departing from a seed represented by a regular 
hexagon. The cells above are obtained by drawing vectors bi of unit length 
and then completing the hexagons such that oposite sides are parallel and of 
the same length. 




Figure 4. Case 1: Two dimensional lattice obtained by setting 9 = 0, \ = 1, 
uj = 1/15. We start by deforming the vector b3 near the seminal cell at the 
origin. For the rest of the lattice, we follow the recurrence relations and the 
construction indicated in the text. Periodicity appears in the direction bl — b2 
indicated with a line at 60 degrees and passing through the origin. 




Figure 5. Case 2: Two dimensional lattice obtained by increasing the angle 
9 = 7r/4. As before, A = 1, cj = 1/15. In this case both vectors b2 and b3 
are deformed near the seminal cell at the origin. The rest of the lattice is 
constructed by using the recurrence relations and completeting the hexagons 
such that opposite sides have equal length. There is no periodic symmetry in 
this case. 




Figure 6. Case 3: Two dimensional lattice obtained by setting 9 = ir/2. As 
before, A = 1, w = 1/15. In this case, the vector b2 is deformed near the 
regular cell at the origin. The rest of the lattice is constructed by using the 
recurrence relations and completeting the hexagons such that opposite sides 
have equal length. There is no periodic symmetry in this case. 
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the cavity. For a cavity of height h (in the z-direction) and for frequencies u < c/2h, 
the scalar field obeys the Helmholtz equation 



For the top and bottom plates, it fulfills Dirichlet or Neumann boundary conditions, for 
TM and TE modes respectively. For cavities with a typical horizontal extension of a tens 
of centimeters and a typical height of a few millimeters, the frequency range of interest 
lies in the domain of microwaves. Such cavities constitute paradigmatic examples for a 
variety of phenomena in open and closed two-dimensional systems as discussed in the 
introduction. 

The experimental set-ups using chaotic microwave cavities adopt various 
configurations depending on the specific studies they are intended for, but the technique 
to feed microwaves into the cavity and to collect the signal of interest is ubiquitous. The 
central device is the network analyzer which performs emission and lock-in detection 
of microwaves from a few tens of MHz to a few tens of GHz. The cavity is linked 
to the network analyzer through flexible coaxial 50i7-cables connected to monopolar 
antennas whose central conductor penetrates into the cavity. Usually, several antennas 
are dispatched over the top and/or bottom plates. For a measurement, only one antenna 
at a time is used as a microwave emitter and another (in transmission) or the same (in 
reflection) as a receiver. The other unused antennas are terminated by 50 fl loads so that 
all antennas behave the same way regarding the losses they imply. The measurements 
are given in terms of scattering coefficients which form the complex ^'-matrix 



where Sn (resp. S22) measures the reflection on port 1 (resp. 2) and 5*12 (resp. 
5*21) measures the transmission from port 2 (resp. 1) to port 1 (resp. 2). All the 
measurements are performed after a proper calibration to get rid of any parasitic 
influence of cables and connectors and even of the analyzer itself. 

In a microwave cavity, a varying potential can be obtained by introducing 
substances of a varying permitivity e(r). The wave equation then reads: 



not preclude the quantum-classical analogy. 

Recently, one of us developed experiments implementing equation (155!) in a 
disordered microwave cavity [37]. The physical phenomenon put under scrutiny was 
Anderson (or strong) localization (see [H] for a recent survey of this prolific domain). A 
network analyzer Rohde & Schwarz ZVA-24 was used in a frequency range from 1 GHz 
to 10 GHz. The disordered potential was introduced through 200 dielectric cylinders 
(Temex-Ceramics, E2000 series) of high dielectric permitivity (e = 37) and low loss 
(quality factor Q = 7000 at 7 GHz). Their height fitted that of the cavity, 5mm, and 



V'V^(r) = e 



(53) 




(54) 
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several diameters ranging from 6 mm to 8 mm were used. The disorder was numerically 
generated and the scatterers were precisely positioned. In this experiment the central 
conductor of the antennas was perpendicular to the plane of the cavity. Thus, the TM 
polarization of the electromagnetic field was selected. [37] [15] . 

The same dielectric cylinders, used for the localization experiments, can be arranged 
in periodic or other ordered patterns. In the domain of optics, periodic structures in 
semiconductor materials are widely used to obtain particular transport properties of 
guided light. Thanks to photonic crystals (also called photonic band gap materials) 
the technology of photons conspire to supplant the one of electrons in domains such 
as communication and information technologies, computing and sensing [9]. Microwave 
cavity experiments, which are definitely intended for more fundamental issues, bear 
the benefits of their versatility. They allow, for example, distortions destroying the 
periodicity of the potential such as the ones described in sections 2 and 3 for emulating 
Dirac oscillators. 

As emphasized in section 3.2, we have to apply the tight binding condition. 
In the localization experiment, the field filled all the cavity, TM polarization being 
supported inside and outside the resonators. For the experiments we proposed with 
ordered structures, the requirements are quite different. In order to transmit the 
energy efficiently into the cylinders, we use in-plane antennas, and consequently TE 
polarization. Due to the wavelength reduction by a factor ^/e ^ 6 inside the dielectric 
material, TE modes can be excited into the scatterers above 5 GHz, while 30 GHz is 
required in air between the resonators. In the proposed frequency range, resonators 
support modes which decay evanescently in the surrounding space. This implies 
exponential decay outside the resonators. For a dielectric cylinder of 8 mm in diameter, 
the first TE mode appears at 6.66 GHz. We experimentally checked that this mode 
is isotropic: ip{r) ~ Jo{r), thus enabling isotropic coupling between resonators. We 
studied the range of the coupling and its spatial dependence. Measurements were 
done with 2 and 3 resonators equally spaced on a line, and 6 resonators placed at the 
vertices of a hexagon (a benzene-like structure). For all configurations, we observed an 
exponentially decreasing coupling with a characteristic length of 400 m~^. The 3-cylinder 
and benzene measurements clearly established that the second-neighbor coupling is 
negligible for distances between the centers of the scatterers above 10 mm. This gives a 
large range of the coupling constant for which the tight-binding model is fulfilled. All 
these experimental results will be published in a forthcoming paper focused on transport 
properties in graphene-like structures [39] . 

We have thus established that we can meet the conditions for the emulation of the 
Dirac oscillator and related problems can be met with arrays of dielectric microwave 
oscillators between two condcuting plates. These conditions will also allow to emulate 
other Dirac operators corresponding to gyroscopes, disordered systems, etc. Note though 
that it will be difficult to emulate a Dirac hydrogen atom as the distances between 
resonators should become extremely small, causing conflict with the diameter of the 
discs. 
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5. Conclusions and outlook 

We have proposed that a wide class of relativistic equations of the Dirac type can 
be implemented by arrays of potential wells. Such quantum systems can be well 
approximated by tight binding hamiltonians if sufficiently fast exponential decay of the 
wave functions in the classically prohibited region is ensured. This idea was especifically 
implemented for the Dirac oscillator. Furthermore we have shown that we can implement 
such a situation experimentally with arrays of dielectric resonators. The use of TE modes 
and the conducting plates below and above of the resonators ensures at appropriate 
frequencies well isolated resonances in the resonators and exponential decay outside. 
We have thus revealed a practical way to emulate Dirac-like equations for both massive 
and massless particles using classical wave systems. 

While we presented an emulation for the Dirac oscillator in one and two dimensions, 
it is clear that our algebraic treatment allows other possibilities. The most promising 
example is the emulation of Dirac gyroscopes |I21J as the number of states in this case 
is finite to begin with due to conservation of total angular momentum. This forces 
its realization on finite grids without approximations in that respect. The relativistic 
hydrogen atom, on the other hand presents obvious difficulties regarding the steep 
potential needed near its singularity and the nearly fiat potential at large distances 
requires separations too small for the radii of the resonators used. Some intermediate 
region of Rydberg states could possibly be emulated. Leaving the realm of integrable 
systems, it might be interesting to introduce random small perturbations in the positions 
of the wells. This might mimic some properties of random matrix Dirac operators jlT] . 

Appendix A. Eigenfunctions 

We first determine the ground state 0o- We write 





m=0 



and use a0o = to obtain the recurrence equation 




(A.2) 



with solution 




(A.3) 



Applying raising operators we obtain 
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(A.4) 



The space-dependent wave functions 0n(r) are obtained by defining a specific form of the 
locahzed functions C,a,C,b- For definitness, we choose these functions to be the ground 
state of a deep square well. The ground state 00(1") and the density |0o(i")P are shown 
in figures lAll and IA2[ 



Ground State 




- X [A] 



Figure Al. Ground state wavefunction in position space (in units of A). The 
localized wave functions are chosen to be the groundstate of individual wells 
with energy a = (3 = 1. The width of the wells is 2A/3. The parameters of the 
lattice are nmax = 20, uj = 1/20 and A = 1. A gaussian envelope is visible. 
The signs alternate from site to site. 




Figure A2. Ground state density as a function of position in units of A. The 
localized wave functions are the groundstates of individual wells with energy 
ct = /3 = 1 and width 2A/3. The parameters of the lattice are JT-rraaz ~ 

20, 

Lo = 1/20 and A = 1. A gaussian envelope is visible. 
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